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I. INTRODUCTION 


Finite state space semi-Markov models find application in a variety of areas such 
as queueing theory, reliability, and clinical trials [Refs. 1,2,3]. The application of these 
models often centers on the distribution of a /first-passage time to a state or a set of 
states representing for example the lifetime of a system or the end of a busy period of a 
server. Suppose that the observations of the path of the semi-Markov process are all 
that is Known about the process. 

In a number of these areas, data arise that are censored. This happens 
frequently, for instance, when fitting lifetime distributions either in medicine or in the 
field of industrial quality control. In medicine, one might be measuring the amount by 
which some new drug extends the life of terminally 1ll patients. A certain number of 
patients are still alive at the end of the experiment, so we do not know how much their 
lives have been extended overall, and certain others might have died of unrelated 
causes or have been removed from treatment prematurely. In quality control one 
might be measuring the distribution of time-to-failure for a sample of integrated circuit 
chips under conditions that accelerate aging. Again, many of the chips may not have 
failed by the end of the trial, while others may have failed at the very beginning due to 
manufacturing defects unrelated to the mechanisms which cause failures in the long 
run. | 

This thesis reports the results of a simulation experiment to compare various 
parametric and nonparametric estimates of the distribution of a first-passage time for a 
particular semi-Markov process with censoring. The specific simulation model and 
estimates considered are given in Chapter 2. Chapter 3 contains the details of the 


simulation experiment and results. Conclusions from the study are given in Chapter 4. 


Il. NATURE OF THE PROBLEM 


A. PROBLEM 

Suppose we observe N individuals. Let X,(i) be the state of the i individual at 
time ¢. We will assume a), t20} i=1, 2, .... N, are independent identicamy 
distributed semi-Markov processes with three states {0,1,2}. The individuals start at 
t=0 in state 1. Upon leaving state I, they transition to state 0 with probability 0 and 
to state 2 with probability 1-8. From state 2, transition is to state 1 with probability 1. 
State 0 is an absorbing state. The sojourn time in state / has a distribution function F. 
(i= 1,2). The individuals are censored, independently. The censoring times are 
exponentially distributed with mean 1/c. The entire path of transitions and sojourn 
times are observed until the time of censoring, if any. Let D be the first entrance time 
to state 0. The problem is to estimate the survival distribution P{D > 2} with the 


censored data of the N individuals. 


B. ESTIMATORS 
Four estimators for P{D > 1} will be described in this section. The first being 
the Kaplan-Meier estimate [Ref. 4], and the others are Maximum Likelihood, Renewal 


Equation, and Asymptotic Renewal estimates from a paper by P. A. Jacobs [Ref. 5]. 


1. Kaplan-Meier Estimate 
One nonparametric estimate for censored data 1s the product limit estimate. 
Let U,, U,, .., U, be independent identically distributed random variables with 
distribution G. Let vr eS ae Ve be independent identically distributed times to 


censorship. Let 


Z, = min (U.,V,) (eqni2) 
and 
OFS el Ot (eqn 2.2) 
1 otherwise . 


ee Z(1) = Z(2) ne Z(n) be the order statistics of {Z.} and 8(;) bemtne 


corresponding order statistic of {6,}. The Kaplan-Meier estimate of G(t) is 


nm ca int < Zin) 


Gt) =< 0 ift > Zn) & 8qy = 0 (eqn 2.3) 
Undefined ites Z(n) & 9(n) = | 
where 
Ca) = (n-i)/(n-1+ 1) (eqn 2.4) 


[Ref. 4:p. 464] and G(t) = [-G(t). If there isn’t any censoring, then the product limit 
reduces to the binomial estimate for each ¢t. This estimate applied to the data of the 
passage times to state 0 for the N individuals will be referred to as the Kaplan-Meier 
estimate of the distribution of the first passage time to state 0 and denoted as 
P(t) = P{D>v. 


2. Maximum Likelihood Estimate 
In this subsection, the maximum likelihood estimate will be given for the 
special case when the sojourn time in state / is exponentially distributed with mean I/p, 
(i= 1,2). 
eet Ri, be the number of transitions from state / to j for one individual. The 


log likelihood function for the individual 1s 
£ = R,,In(1-9)+ R, gin 9+ R,,1n p,+(Rjg + Rj) In p, - P,T, - PT, (eqn 2.5) 


where T, (i=1,2) is the total time spent in state / before entrance into state 0 or 


censoring [Ref. 5:p. 2]. The maximum likelihood estimators are 


A 
9 = Rig (Rig t Ry >) (eqn 2.6) 


Pp, = (RyptRpo)/T, (eqinez 7) 
A 
prea ae lo (eqn 2.8) 


The maximum likelihood estimate for the survival distribution 1s 
(Ren ope ceqnul. 17] 


A An A A A A A A A A _ A A 
P {D> t} = {0p /(A,-Aa)} [Ay + P9)/Ag|explta,]-[(A, + py)/A,Jexp[taj]} (eqn 2.9) 
A AN 
where 4, and A, are the roots of the equation 
AAA A A 5 
Os ae (Bip Day) teeta (eqn 2:10) 


The above estimate will be referred to as the parametric estimate and denoted as 
A A 
Pt) = PAiD> vw. 


3. Renewal Equation Estimate 


The probability P{D > tf} satisfies the renewal equation 
P{D>t} =F (t)+(1-0)f)F, (ds)F,(t-s) + (1-6) (EF *F ds)PiDe ts} (eqn 23) 
where F, is the distribution of the sojourn time in state i, F, (t) = I-F, (t), and Fine 


the convolution of er and F.. 


The solution to the renewal equation 2.11 1s 


P{D>t} = g(t) + |’, R(ds) g(t-s) (eqn 2.12) 
where 
g(t) = F(t) + (1-0) fy F)(ds) F,(t-s) (eqn 2.13) 


10 


and 


RO) = ¥ 1-6)" (FSF) (eqn 2.14) 
mail 
where (F,*F,)" (t) denotes the n-fold convolution of (F,*F.,) with itself at time ¢. 

A nonparametric estimate for P{D > 1} can be obtained by replacing F. by its 
Kaplan-Meier estimate and @ by its maximum likelihood estimate in equation 2.12. If 
the largest sojourn time in state / 1s censored then the Kaplan-Meier estimate of F. is 
not an honest distribution function (F,(00) <1) since the estimate is undefined past the 
largest sojourn time. In this case the dishonest distribution estimate 1s used in all the 
remaining computations which will give a conservative estimate of the survival 
distribution. 

An approximation to equation 2.12 can be found by using a discrete time 


approximation to R(t) as follows. Let 6 > 0 be a constant and let 
p,(0) = (1-8) {[F,*F,](nd) - [F,*F,]((n-1)0)} . (eqn 2.15) 
Recursively approximate R(t) as follows 


in) = 0 (eqn 2.16) 
R,(0) a p, (6) 


R,(n6) = » p,.(5) + p,. (8) R,((n-k)6) . 
=f 


An approximation to the solution of equation 2.12 using estimates of F. and 0 1s 


mle) 
P (D> t} = g(t) + 5 (R,(kd) - R,((k-1)8)} 2(t-Kd) (eqn 2.17) 
R= 
where n(6) is the largest integer less than t/6 [Ref. 5:p. 9 eqn 2.9]. If the number of 
individuals N or the time ¢ are large, the estimate of equation 2.17 may require a large 
number of additions of small non-negative numbers. This estimate will be referred to 


“A 
as the renewal estimate and denoted as P (t) some 
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4. Asymptotic Renewal Estimate 
A A 
Let F. be the Kaplan-Meier estimate of F, and 0 be the maximum likelihood 


estimate of @:; then define 


A A 
(5) = 5 expls&] F.(ds) (eqn 2.18) 
A 
where again F. may be a dishonest distribution due to censoring of the last sojourn 
time in state i. The asymptotic renewal estimate of the survival distribution 1s 
[Retss: 0.81 seqes3- i 


P (D>t} = expltt] (b/f) (eqn 2.19) 


A ee 
where K 1s the solution to the equation 


A A A_A 

(1-8) (K) $,(K) = 1 (eqn 2.20) 
and 

A A A A fn | 

uw = (1-0) fo explsk] s (F, *F,)(ds) (eqn 2.21) 
and 

“™ AA A OA 

b = (6/K) @,(K) . (eqn 2.22) 


A 
The « for equation 2.19 was found by numerical search using equations 2.18 and 2.20. 
The above estimate will be referred to as the asymptotic estimate and denoted as 
A A 
a) eee, 
“w 

If P{D>t} were exactly the solution of the equation 2.12 with the 
Kaplan-Meier estimate of F. and the maximum likelihood estimate of 8 being used 
then 


A N 
PAD i Daal (eqn 2.23) 


as [-»° in the case where the Kaplan-Meier estimates are honest distributions. 


II. ANALYSIS OF THE PROBLEM 


A. SIMULATION 

A Fortran program is written to generate and analyze the data for this problem. 
All simulations are carried out on an IBM 3033AP computer at the Naval 
Postgraduate School using the LLRANDOM II random number generating package 
[Ref. 6]. The data for the simulation experiments are generated as follows: 
Independent exponential censor times with mean I/c are generated for each individual. 
The individual starts in state 1 at t=0 and an exponential time with mean 1/p, 1s 
generated for the sojourn time. A comparison between the sojourn and censor time is 
done; if the sojourn time is smaller, then the sojourn time 1s recorded; if the censor time 
is smaller, the truncated sojourn time and the censored death time are recorded. From 
state 1, if not censored yet, a uniform random number 1s compared to theta; if less 
than theta, the process jumps to state 0 and the uncensored death time is recorded; if 
greater than theta, the process jumps to state 2 and an exponential sojourn time with 
mean 1/p., is computed. The total time (sojourn times in state 1 plus sojourn times in 
State 2) is compared to the censored time; with the same actions as listed above. From 
state 2, the process jumps to state 1 and continues until an uncensored or censored 
death occurs. The times are recorded and the next individual 1s started. This continues 
until all N individuals have been generated. The data in each state 1s sorted in 
increasing order for ease of program manipulations. If N is small, it is possible for all 
the sojourn times in a state to be censored or for all the first passage times to state 0 to 
be censored which results in P(t), P(t), or P, (t) being undefined for all ¢. In these 
cases the replication is dropped and a new replication generated. 

A sample data set is listed below for N= 10. The first row under state | and state 
2 gives each particular censored or uncensored sojourn time that is generated for that 
state. Under each sojourn time, the binary number indicates whether the individual is 
censored (1) or not (0) during that sojourn time. State 0 indicates times of death 
(passage time to state 0), and whether censored (1) or not (0); note that the times 
indicate either the time of death (not censored) or the time of censoring (censored 
death time). The sojourn and death times listed below have been sorted, along with its 


associated censor indicator. 


13 


State 


0.1629 0.2041 2.2201 
0 0 0 


State | 


0.1356 0.1615 0.2114 0.2748 0.2996 0.3067 0.3450 0.3725 0.3996 0.4305 


l ] ] 0 0 0 I 0 0 I 
0.8676 1.1980 2.4630 
0 0 0 


StateO N=10 


0.1356 0.1615 0.2748 0.3450 0.8676 0.8832 0.9930 1.1980 2.4630 2.7312 
] l 0 I 0 0 ] 0 0 I 


Using equations 2.3, 2.9, 2.17, and 2.19, estimates of the survival distribution 
P{D>t} from the data are calculated from subroutines in the Fortran program. Output 
from the program produces a table like the one below that includes: time, actual 
suv val probability (ACT(t)), ee estimate (P, {D>t}), renewal estimate 
(P. {D> t}), asymptotic estimate (P, {D> t}), and the Kaplaatieier estimate of the first 
passage time to State 0 8, (D>). The actual survival probability ACT(t) is 
computed using equations 2.9 and 2.10 with the actual parameter values instead of the 
estimated values. The Kaplan-Meier estimate uses only the uncensored first passage 
times to Sent: Output in Table I is for the data set listed above. 

In Table I, the renewal and asymptotic estimates decrease as ¢ increases. In this 
case, the largest sojourn times in both state 1 and state 2 are uncensored. To 
demonstrate what can happen when the largest sojourn times are censored, Table I] 
shows a case where the largest sojourn times in state I and state 2 are censored. 
Notice that after t=5 there is little change in the renewal estimate. The survival 
probability levels off and becomes constant. The asymptotic estimate starts low (half 
the probability) and goes to zero just after t=5. In a third case, when either of the 


largest sojourn times in state 1 or state 2 are censored, the effects are somewhere 
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between the two cases mentioned above; the renewal estimate starts to level off but 
may not become constant and the asymptotic estimate starts lower than normal and 
may go to zero. The dishonest Kaplan-Meier estimate of F. has a definite affect on 
A A 

P(t) and P(t) for large ¢. 


AB Ube 
OUTPUT FROM PROGRAM 


Survival Probability P{D> t} 








Time _ ACT(t) P (D>) P{D>t}  P,{D>t} P(D> 4) 

= 0.79965 0.73641 0.66667 0.68473 0.87500 

1.0 0.66340 56522 0.52606 0.54348 0.58333 

2.0 0.47996 0.35318 0.36158 0.34238 0.38889 
0 0.19737 0.09549 0.09181 0.08560 Undefined 
7.0 0.10985 0.04027 0.03451 0.03397 Undefined 
10.0 0.04563 0.01103 0.00874 0.00849 Undefined 
ID 0.02194 0.00375 0.00272 0.00268 Undefined 
15.0 0.01055 0.00128 0.00086 0.00084 Undefined 

TABLE II 
OUTPUT FROM PROGRAM 
Survival Probability P{D >t} (largest sojourn censored) 
A wN “A ” 

Time ACT(t) Pet} PDEA Payor; Pe ty 

5 0.79965 0.79130 0.78783 Gi35520 0.90000 

1.0 0.66340 0.65115 0.65547 0.16873 0.78750 
2.0 0.47996 0.46358 0.52706 0.03851 Undefined 
5.0 0.19737 0.18049 0.49071 0.00046 Undefined 
7.0 0.10985 0.09675 0.49020 0.00002 Undefined 
10.0 0.04563 0.03797 0.49017 0.00000 Undefined 
IS 0.02194 0.01742 0.49017 0.00000 Undefined 
15.0 0.01055 0.00799 0.49017 0.00000 Undefined 


15 


B. ANALYSIS 

For the simulated model described above, parameter values of p= ie p= i 
0=0.5, and c=0.5 are used. The simulation uses two different numbers of observed 
individuals. The number of individuals is set at 10 and 50, representing a low and 
moderate number of subjects. The simulation is replicated 500 times utilizing different 


seeds to generate the data. The average relative bias for each estimate is computed by 


ARB(t) = (1/M) s (EST,(t)-ACT(t))/ACT(t) (eqn aay 
A*f 
where EST.(t) is the value of an estimate computed for the i replication at time ¢ and 
ACT(t) is the actual model value at time ¢. For the Kaplan-Meier estimate, M is taken 
as the number of Kaplan-Meier estimates of the first passage time to state 0 still 
defined by time ¢. For the other estimates, M 1s the number of replications (500). 

The figures below show histograms of the relative bias of the observations 
(EST (t)-ACT(t))/ACT(t). Figure 3.1a shows histograms of the relative bias for each of 
the four estimates when N= 10 and at r=0.5. Each of the histograms looks relatively 
normal with possibly a slight skew to the left. The parametric estimate has the tightest 
distribution and the asymptotic estimate the worst which is expected since the 
asymptotic properties are for large ¢. Figure 3.1b shows the relative bias for each 
estimate when N=10 and r=5.0. The parametric is somewhat normal but skewed to 
the right. The renewal estimate looks a little less skewed. The asymptotic estimate is 
skewed to the right and looks exponential. At time r=5.0, less than half of the 
Kaplan-Meier estimate of the first passage time to state 0 are defined. The histogram 
of the defined Kaplan-Mcier estimate is starting to show an accumulation of mass at 
-1.0 which is the value of the relative bias where the largest passage time observation 1s 


uncensored and less than 5.0. 
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Figure 3.1a Histograms of relative bias for N= 10 and t=0.5. 
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Figure 3.1b Histograms of relative bias for N= 10 and t=5.0. 
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Table III shows the ARB(t) of the estimates for the case when N= 10 individuals 
and Table IV for the case when N= SO individuals. The ARB(t) for each estimate is 
given for selected values of ¢. Along with the ARB(t) in the parentheses is the 
corresponding standard error. The standard error is computed by taking each 
observation of the relative bias (EST.(t)-ACT(t))/ACT(t) and subtracting the ARB(t), 
Squaring this and summing over all M observations, then dividing by M-1l. This 
produces the distribution variance, which is divided by M and the square root taken of 
to get the standard error of the ARB(t) for each estimate at time ¢. The variance 
together with the average relative bias can be used to obtain an estimate of the relative 
mean squared error of the estimate. The right most column of the Tables III and IV 
gives the number of replications out of 500 that still has defined Kaplan-Meier 


estimates of the distribution of the first passage time to state 0 by time rf. 


TAR EI 
BVEKAGE REEATIVE BIAS 


Exponential Model N=10 (500 Reps) 

















: IS A A A 
Time P(t) et) P (t) Pg) HK M 
p r a k 

a, -.00183 01292 -. 18069 0144 500 
(00533) (00545) (01287) (00743) 

1.0 00575 02788 -.24628 04135 499 
(00942) (00917) (01501) (01153) 

2.0 .03470 .05672 -.34206 .07631 462 
(01656) (01568) (01942) (01847) 

5:0 4324 mS Gn -.43114 -.37581 225 
(.04126) (04267) (03196) (07156) 

7.0 047 1.30581 -.39861 -.69076 185 
(.06596) (07867) (:04399) (08469) 

10.0 1.15728 3.84602 -.25101 - 1.00000 174 
(12758) (.19504) (07358) (00000) 

WS 2.06816 8.49575 -.02144 - 1.00000 174 
(.21894) (.41214) (.11647) (00000) 

15.0 3.52449 8.05457 .34887 - 1.00000 174 
(37592) (.86673) (18841) (00000) 


Ne 


The parametric estimate P(t) uses the most correct information about the 
process. For N=10, the parametric estimate is within three standard deviations of 
zero bias forr<5. As ¢ gets larger, the relative bias tends to increase. The parametric 
estimate understandably has the smallest relative bias for small ¢. For large r, the small 
sample sizes involved are probably responsible for the larger relative bias. For small 
times the renewal estimate and the Kaplan-Meier estimate for the distribution of the 
first passage time to state 0 have about the same average relative bias. For sinall N 
and large r, the renewal estimate has large bias. As noted before, the renewal estimate 
will be biased if the largest observations of the sojourn times in a state are censored 
thus causing the Kaplan-Meier estimate F, to be undefined. The bias could also be 
caused by the step size in the discrete time approximation (step size 0.01) being too 
large, or by numerical error in summing large quantities of small numbers, as 
mentioned earlier. The Kaplan-Meier estimate does well for small t and small N. As 
time increases, the number of data points depreciates rapidly. Because of the small 
number of subjects in each run, the Kaplan-Meier estimate of the distribution of the 
first passage time to state O lost over half its data due to undefined distributions. By 
time r= /0, there are no survivors using the Kaplan-Meier estimate, resulting in the 
-1.0 average relative bias. From equation 2.23, the renewal estimate and the 
asymptotic estimate should be approximately the same for large rif the Kaplan-Meier 
estimates F, are always defined. The asymptotic estimate is negatively biased for small 
t but changed over at t> /2.5. Once again, it could be biased due to censoring of the 
largest sojourn times. The asymptotic estimate has the smallest average relative bias 
for large time ¢. 

Figure 3.2a snows histograms of the relative bias for each of the four estimates 
when N= 50 and at r=0.5. Each of the histograms again looks relatively normal.The 
distributions are much tighter when compared to Figure 3.la. The parametric estimate 
has the tightest distribution and again the asymptotic estimate the worst. Figure 3.2b 
shows the relative bias for cach estimate when N=50 and :=5.0. All the estimates 
except the Kaplan-Meier estimate look relatively normal with possibly a slight right 
skew. The Kaplan-Meier estimate of the first passage time to state 0 has just over two 
thirds of its distributions defined and is showing the start of an accumulation at -1.0 


due to the largest passage time to state 0 being less than 5.0. 
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50 and t=0.5. 


Figure 3.2a Histograms of relative bias for N 
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50 and t= 5.0. 


Figure 3.2b Histograms of relative bias for N 
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TABLE IV 
AVERAGE RELATIVE BIAS 


Exponential Model N=50 (500 Reps) 




















A A A A 2 
Time P pit) P(t) a) P(t) #KM 
LD -.00430 -.02928 -. 12087 -.00082 500 
(.00398) (.00455) (.00743) (.00512) 
1.0 -.00567 -.11500 -.12875 00403 500 
(.00220) (.00278) (.00580) (.00330) 
2.0 -.00364 -.28952 -. 17511 -.00003 500 
(.00705) (.00742) (.01087) (.00872) 
0 03428 -.03348 -.29824 01067 34] 
(.01623) (.01946) (.01876) (.04022) 
70 08483 Goes -.33861 -.35364 224 
(.02326) (.03139) (.02347) (.07749) 
10.0 oe. 2.01703 -.35577 -.75760 188 
(.03633) (.07707) (.03175) (.08671) 
PES Some) 4.36834 -.33896 -.72166 186 
(.05091) (.16959) (.04109) (.14855) 
15.0 52199 9.2217] -.29647 - 1.00000 184 
(.07067) (2307 (.05409) (.00000) 


For N=50, the parametric estimate again does well for small ¢. The average 
relative bias is within three standard deviations of zero for <7. An improvement for 
large ¢ for the parametric estimate is expected because of the increased number of 
individuals. The renewal estimate again has large bias for large ¢, though not as much. 
The asymptotic estimate is negatively biased throughout time. For ¢> 2, the bias looks 
constant. For N=50, the Kaplan-Meier estimate of the first passage time to state 0 
has negligible average relative bias for t< 7. However as ¢ increases, the Kaplan-Meier 
estimate looses an appreciable amount of its data due to undefined distributions. By 
t= 15, the Kaplan-Meier estimate has no survivors. Once again for large ¢ (t= /5), the 


asymptotic estimate has the smallest average relative bias. 


os 


A simulation experiment is done for a case in which there is a relatively high 
number of individuals N= 100; c=0.5, p,)=1, p,=1, 9=0.5. The results appear in 
* Table V. The increased number of individuals has decreased the average relative bias 
for all the estimates. The standard error of all the estimates has also decreased. From 
Tables III, 1V, and V, it appears that as the number of observed individuals increase 


the average relative bias for all the estimates decrease. 


TABLE V 
AVERAGE RELATIVE BIAS 


Exponential Model N=100 (500 Reps) 




















A “A “A . 
Time E a t) Pv Paw P(t) #K M 

-.00192 -.09574 -.09600 -.00379 500 
(.00150) (.00251) (.00399) (.00240) 

1.0 -.0027 1 -.13534 -.08705 -.00235 500 
(.00275) (.00318) (.00512) (.00367) 

2.0 -.00256 -.38676 -. 10980 -.00267 500 
(.00495) (.00426) (00756) (.00609) 

5.0 01258 -.69527 -.21245 04573 448 
(.01136) (.00808) (.01378) (.02474) 

7.0 03523 -.51234 -.26003 = Soa 281 
(.01588) (.01557) (.01729) CO>v) 

10.0 08765 28800 -.30374 -.47180 200 
(.02337) (04112) (.02242) (.11283) 

he) Ads 1.96444 -.32021 -.87080 185 
(.03059) (.09357) (.02708) (.09137) 

15.0 22930 5.69514 -.32195 -.87581 185 
(.03909) (.19744) (.03245) (.12419) 


In order to investigate the effect of censoring on the values of P(t) and P (t) for 
large ft, a simulation study is done in which the exponential censoring tunes has a mean 
of 1/c= 1000. The other parameters are p,=1, p,=1, and 0=0.5 as before. Once 
again the number of individuals are N= 10, 50. The results are presented in Tables VI 


and VII. The results in Table VI suggest that the effect of the small sample size 
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resulting from N=10 dominates the performance of all the estimates except for the 
Kaplan-Meier estimate for large z. The results of Table VII suggest that the method of 
computing P -(t) is affecting its performance for large ¢ since P (w/P Bee Was => 00! 
For N= 50, limited censoring has improved the average relative bias of all the estimates 
for large ¢«. Somewhat surprisingly, with limited censoring the Kaplan-Meier estimate 
has almost the best average relative bias. However, the standard error of its estimate is 
larger than that of the other estimates’. Thus the Kaplan-Meier estimate tends to be 


more variable than the other estimates. 


TABLE VI 
AVERAGE RELATIVE BIAS 


Exponential Model N=10 (Limited Censoring, c=0.001) (500 Reps) 

















: A “A A “A . 
Time P(t) P (t) P(t) P(t) #KM 
-.01895 -.01133 -.03853 -.00768 500 
(:00360) (:00441) (:00476) (:00718) | 
1.0 -.03008 -.04266 -.02281 -.00492 500 
(:00629) (:00644) (00614) (:01025) 
2.0 -.04178 -.08352 -.02283 -. O1816 500 
(:01058) (:00955) (:01049) (01461) 
5.0 -.01755 -.02716 -.01934 04571 500 
(:02166) (:02091) (:02207) (:02963) 
7.0 .04003 05500 02483 09415 500 
(:02968) (03044) (:03009) (03971) 
10.0 .18364 Pe .15076 05896 500 
(.04466) (.04946) (.0448 1) (.0642 1) 
> .36140 33159 Suey 03902 500 
(.06146) (.07442) (.061 27) (6093535) 
15.0 .60294 94267 ou oe -.00407 500 
(:0841 4) (11360) (08351) (113365) 
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TABLE yal 
AVERAGE RELAFIVE BIAS 


Exponential Model N=50 (Limited Censoring, c=0.001) (500 Reps) 




















A A A A 
Time Pt) Ee eae) P,(t) #KM 

-.00470 -.04348 -.07253 00380 500 
(.00144) (.00193) (.00267) (.003 13) 

1.0 -.00760 -. 15247 -.03784 00104 500 
(.00259) (.00285) (.00318) (.00425) 

2.0 -.01110 -.40891 -.01343 00562 500 
(.00452) (.00355) (.00464) (.00645) 

5.0 -.01035 -.79859 -.00883 00497 500 
(.00982) (.00278) (.00987) (.01279) 

7.0 -.00123 -.79946 -.00283 -.003 50 500 
(.01340) (.00241) (.01350) (.01820) 

10.0 02459 -.59136 01905 -.01116 500 
(.01903) (.00335) (.01922) (.02926) 

[2 05757 -.13724 04945 -.08114 500 
(.02413) (.00679) (.02440) (.04066) 

15.0 NOMS 84405 09097 -.06882 500 
(.02977) (.01457) (.03015) (.0563 1) 


Below are reported simulation results experimenting with different parameter 
values of p, and c. For these studies, the number of individuals is set at N=50 to 
reduce the effects of undefined Kaplan-Meier estimates of F.. Four different cases are 
simulated. The sojourn time in state 1 is changed to reflect a higher and lower mean 
sojourn time and the censoring mean time is changed to reflect more or less censoring. 

The first cases that are simulated are the changes in the mean sojourn time in 
state 1. The other parameters are c=0.5, p,=1, and 0=0.5 as before. The mean 
sojourn time of state 1 is increased from 1 to 2 (p,=0.5) and decreased from | to 0.5 
(p,=2). With the increase in the mean sojourn time of state 1, the probability of a 
death being censored increases. For a decrease in the mean sojourn time, the opposite 
is true. There are quicker jumps out of state 1, resulting in more uncensored deaths. 
Tables VIII and IX show the computed average rclative bias using equation 3.1 along 


with the associated standard error. 
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TABLE VIII 
AVERAGE RELATIVE BIAS 


Exponential Model N=50 (p,=0.5) (500 Reps) 




















; A AN A A 
Time P(t) P(t) P(t) P(t) #KM 

oD 00045 -.QOSO5 -. 14754 .00062 500 
(.00136) (.00188) (.00921) (.00234) 

1.0 OO175 -. 02090 -.18477 .00707 - 500 
(.00261) (.00326) (.01020) (.00366) 

2.0 064 -.06266 -.29529 00877 500 
(.00496) (.00597) (.01228) (.00616) 

5.0 03886 .06045 -.40236 1g 415 
(.01198) (.01328) (01702) (.02202) 

0 07539 .29430 -.45897 -.04033 Zoo 
(.01717) (.02093) (.01959) (.05350) 

10.0 15280 0927 -.50727 -.65507 ey 
(.02643) (.04165) (.02366) (.07690) 

2-5 .24079 1.53886 -.52490 -.88639 137 
(.03615) (on) (.02772) (.06656) 

15.0 (Soo) 2.19204 s027358 -.86033 7 
(.04857) (.12186) (.03276) (.09855) 


In Table VIII where the mean sojourn time in state 1 increases, the average 
relative bias of the parametric estimate looks about the same as in Table IV. The 
average relative bias of the renewal estimate is slightly better than in Table IV. The 
average relative bias of the asymptotic estimate looks like it increased, but is within 
three standard errors of Table IV. The average relative bias of the Kaplan-Meier 
estimate of the first passage time to state 0 looks the same as in Table IV. The number 
of defined Kaplan-Meier estimates has decreased due to the increase in the probability 
of a censored death as mentioned earlier. There are two survivors at t= /5. 

In Table IX where the mean sojourn time in state 1 decreases, the parametric, 
asymptotic, and Kaplan-Meier estimates have the same average relative bias as in 
Table IV. The number of defined Kaplan-Meier estimates has increased due to the 
decrease in the probability of a censored death. The renewal estimate has increased as 


compared to Table IV. 
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TABLE IX 
AVERAGE RELATIVE BIAS 


Exponential Model N=50 (p,=2) (500 Reps) 




















“\ A A A 
lame P(t) faa) 2 a) P_(t) #KM 
5 -.0069 1 Tots -.11848 00049 500 
(.00323) (.00369) (.00504) (.00434) 
1.0 -.00907 22919 -, 10384 . -.00438 500 
(.00537) (.00512) (.00753) (.00685) 
2.0 -.00826 -.49325 -.16288 -.00674 500 
(.00893) {.00721) (.01187) (O12ZIN) 
5.0 04037 -. 16543 -.30815 -.08353 308 
(.02060) CO2TGy) (.02132) (.05516) 
7.0 11149 95936 -. 34515 -.61487 221 
(.03055) (.04897) (.02811) (.08054) 
10.0 .28446 5.05635 -.34050 -./8435 207 
(.05190) (.15416) (.04336) (161) 
12.5 0557 [2.93165 -.28668 - 1.00000 207 
(.07987) (.40802) (.06490) (.00000) 
15.0. 82543 31.25238 -.17898 - 1.00000 200) 
(12395) (1.11064) (. 10098) (.00000) 


The next cases that are simulated are the changes in the censoring distribution; 
the other parameters are p)=1, p,=1, and §=0.5 as before. The exponential mean 
time to censor is increased from 2 to 4 (c=0.25) and decreased from 2 to | (c=1). 
With an increase in the mean censoring time, the probability of a censored death 
decreases. With a decrease in the mean censoring time, the opposite is true. Tables X 
and XI show the average relative bias for each simulation along with the standard 
error. 

Table X where the mean censoring time increases (1/c=4), falls between Table 
IV and Table VII. The average relative bias of the parametric estimate is worse than in 
the limited censoring case of Table VII but slightly better than in the case c=0.5 of 
Table 1V. The average relative bias of the renewal estimate is much worse than it 1s 


with the limited censoring of Table VII but about the same to slightly better in the tail 
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than in the case c=0.5 of Table IV. The average relative bias of the asymptotic 
estimate is about the same as the limited censoring case but much better than in the 
case c=0.5. The average relative bias of the Kaplan-Meier estimate of the first passage 
time to state 0 is about the same for small to moderate times and worse for large times 
than in the limited censoring case of Table VII and better for large times than in the 
case c=0.5 of Table IV. The number of defined Kaplan-Meier estimates is between the 


two tables, due to the increase in the mean censoring times. There are five survivors 




















past t= 1), 
TABLE X 
AVERAGE RELATIVE BIAS 
Exponential Model N=50 (c=0.25) (500 Reps) 
A A A A 7 
Time P(t) Fa) Pw) P(t) #KM 

= -.00099 -.03349 -.06677 .00 170 500 
(.00179) (.00230) (.00331) (.00325) 

1.0 -.00071 -.13166 -.04538 00874 500 
(.00326) (.00348) (.00434) (.00466) 

ZU 00178 -.36077 -.04544 00837 500 
(.0058 1) (.00525) (.00694) (.00720) 

5.0 02794 -.57486 -.09177 01752 486 
(.01330) (.00892) (.01460) (.02043) 

7.0 .06206 -.27543 -. 10565 06881 4\1 
(.01876) (.01723) (.01949) (.04235) 

10.0 13879 85619 -.09864 -. 10377 307 
(.02822) COs S72) (02727) (.09322) 

5 .2288 | | 2.71999 -.06942 -.59309 276 
(.03788) (.07635) (.03479) (.11491) 

15:0 34587 5.92070 -.01964 -.69629 21D 
(.04987) (.15972) (.04386) (.16349) 


Ze 


In Table XI where the mean censoring time decreases, the average relative bias of 
the parametric estimate 1s about the same for small to moderate times and then is 
worse for large times than in Table IV. The average relative bias of the renewal and 
asymptotic estimates are both worse for r> 2 due possibly to an increase in the number 
of dishonest Kaplan-Meier estimates of F.. The average relative bias of the the 
Kaplan-Meier estimate of the first passage time to state 0 is worse for t>5. The 
number of defined Kaplan-Meier estimates has decreased reflecting the decrease in the 


mean time to censoring. 


TABLE XI 
AVERAGE KEVAI Ears 


Exponential Model N=50 (c=1) (500 Reps) 




















“A “A A a 
Time P(t) P (t) P (t) Pay #KM 
p r a 

.S -.00268 -.01350 = 21295 -.Q00012 500 
(.00270) (.00317) (.00971) (.00369) 

EO -. 00193 -.05220 -.35347 00391 500 
(.00491) (.00568) (.01171) (.00610) 

2.0 00809 -.07041 -.47848 -.00203 488 
(.00878) (.01133) (.01452) COisa2) 

500 .10233 4674 -.66059 -.63062 128 
(.02127) (.03027) (.01795) (07727) 

70 .21600 1.46539 -./0866 =, 93600 joe 
(.03220) (.06044) (.02021) (.05174) 

10.0 47564 4.48047 -./3703 - 1.00000 L11 
(.05588) (.15766) (.02540) (.00000) 

25 SO113 10.05fo2 “(30 Ul - 1.00000 111 
(.08693) (.34001) (03213) (.00000) 

15.0 l20 239 21.62486 -./2108 - 1.00000 111 
(.13558) (.72038) (04211) (.00000) 
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C. ROBUSTNESS 

In the above simulations the maximum likelihood estimate used the known 
correct model. Often, a model needs to be chosen to describe a data set. Attempts are 
made to analyze the data to determine a good model. However, when sample sizes are 
small, the difficulty of finding a good model increases. Hence due to small sample sizes 
or ease of computation, an incorrect model may be chosen to describe a data set. In 
this section, the robustness of the estimates proposed in Chapter II is studied with 
respect to an incorrect model assumption concerning the sojourn time in state 1. 

The data for the simulation experiment in this section are generated from the 
following three state semi-Markov process: Individuals start in state 1 at t=0. The 
probability of a jump to state 0 is 9; to state 2 is 1-@. From state 2 the probability ofa 
jump to state I is 1. State 0 is an absorbing state. The sojourn time in state 2 is 
exponential with mean 1/p,. The sojourn time in state I is the sum of two independent 
exponentials with means 1/p, and 1/p,; that’ 1s, the sojourn time in state | has a 
hypoexponential distribution. Censoring is independent and exponentially distributed 
with mean I/c. The same basic Fortran program is employed, modified for the above 
change. The data generated are analyzed by the same Fortran subroutines for each 
estimate as in the first section. In particular, the (incorrect) maximum likelihood 
estimate of equation 2.9 is used. This maximum likelihood estimate assumes the 
sojourn time in state 1 has an exponential distribution rather than the true 
hypoexponential distribution. 

For the first simulation results reported, parameter values of p) = 1, p,=1, p,;=1, 
§=0.5, and c=0.5 are used. Again, two different numbers of observed individuals are 
used, 10 and 50. The simulation is replicated 500 times and the average relative bias is 
computed utilizing equation 3.1. For the Kaplan-Meier estimate, M is taken as the 
number of defined Kaplan-Meier estimates of the first passage time to state 0 by time 
t. For the others, M is the number of replications. The actual value of the survivor 
function is computed by inverting the Laplace transform of the passage time to state 0 


for the semi-Markov process. 
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Tables XII and XIII show the average relative bias of the hypoexponential 
model at selected values of t along with its associated standard error for N= 10 and S50. 
Again the right most column is the number of defined Kaplan-Meier estimates of the 


first passage time to state 0 out of 500 replications. 


TABLE XII 
AVERAGE RELATIVE BIAS 


Hypoexponential Model N=10 (500 Reps) 




















; A N A A : 
Time E A t) P (t) Pa P, (t) #KM 

e -.04965 00603 -.18749 00434 500 
(.00215) (.00225) (01511) (.00332) 

1.0 -.04224 00518 2275S 00119 499 
(.00415) (.00462) (.01478) (.00622) 

20 02344 1020 -.35821 03759 488 
(.00829) (.00994) (.01670) (.01170) 

5.0 27429 23906 -.45864 -.03008 255 
(.02259) (.02654) (.02448) (.05463) 

7.0 .49966 58490 -.46202 -.51308 164 
(.03501) (.04261) (.03062) (.08072) 

10.0 .98106 hooto5 -.40812 -.90325 142 
(.06148) (.08272) (.04310) (.05049) 

12.5 1.56115 2.93742 <1 120 -.97945 140 
(.09466) (.14238) (.05800) (.02055) 

15.0 2.37288 5.27240 -.18077 - 1.00000 139 
(.14350) (.24293) (.07897) (.00000) 


In Table XII, for N=10, the parametric estimate based on the incorrect model 
shows more relative bias for small ¢ than the results in Table III using the correct 
maximum likelihood model. However, for moderate times /<t</7 the effect of the 
small number of individuals has overwhelmed the effect of the incorrect model and the 


relative bias is approximately the same as for the correct model given in Table III. 
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The average relative bias of the nonparametric estimates appear to do well. The 
renewal estimate and the Kaplan-Meier estimate of the first passage time to state 0 
seem to do very well for small times and about the same for moderate to large times, 
with the Kaplan-Meier decreasing to -1.0 at t= 15. The asymptotic estimate seems to 
do about the same as in the situation of Table III; it is still negatively biased and has 


the smallest average relative bias for large times. 


TABLE XII] 
AVERAGE RELATIVE BIAS 


Hypoexponential Model N=50 (500 Reps) 




















A A A A 
Time P(t) P(t) P(t) P(t) #KM 

5 -.04769 00196 -.10608 -.00089 500 
(:00104) (00104) (00701) (:00154) 

1.0 -.04010 -.00130 -.14888 001 500 
(:00203) (:00229) (00745) (00297) 

2.0 02140 -.02910 -.19231 00557 500 
(100415) (:00520) (:00965) (:00563) 

5.0 21701 05244 -.29425 03315 421 
(01169) (01291) (:01620) (:02328) 

7.0 36307 20321 -.33710 -.09714 246 
(01803) (01998) (01987) (:05655) 

10.0 63698 64409 -.36628 -.53280 162 
(103053) (103839) (02550) ('09457) 

12.5 92760 1.32451 -.36722 -.93675 150 
(:04476) (:06729) (03096) (0524?) 

15.0 1.29131 2.47259 -.35081 -.97947 149 
(06391) (11708) (03768) (:02053) 


In Table XIII, the case of the larger number of individuals N= 50, the effect of 
the incorrect model of the maximum likelihood estimate has a more noticable effect on 
the average relative bias; the average relative bias for the parametric estimate is 
significantly higher than for the nonparametric renewal and Kaplan-Meier estimates for 
t=2. The nonparametric renewal and Kaplan-Meier estimates have about the same 
average relative bias for tS5. The average relative bias of the asymptotic estimate 
does the same as in Table IV, and still consistently negatively biased. The 
Kaplan-Meier estimate of the first passage time to state 0 does the same as in Table 


IV. There is still one survivor at t= /5. 
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A simulation experiment is done for the case in which there is a relatively high 
number of subjects N= 100; the other parameters are as before. Table XIV show the 
effects of the increase in observed individuals. The average relative bias of the 
parametric estimate shows less relative bias than in Table XIII but still significantly 
higher than Table V using the correct model with comparable number of subjects. The 
average relative bias of the nonparametric estimates has lower relative bias than Table 
XIII and about the same relative bias as Table V. Again it appears that the average 


relative bias for all the estimates decrease as the number of individuals increase. 


TABLE XIV 
AVERAGE RELATIVE BIAS 


Hypoexponential Model N= 100 (500 Reps) 




















A A A A E 
Time P(t) P(t) P(t) P_(t) #KM 
5 .00010 04879 -.04324 05034 500 
(.02439) (.02528) (.01982) (.02493) 
1.0 03696 05200 -.05118 07747 500 
(.03963) (.03902) (.03019) (.03903) 
2.0 13258 -.05690 -.05490 .11607 500 
(.05843) (.04927) (.03941) (.05859) 
0 .18619 -.4492] -.21495 -.02024 491 
(.00875) (.01064) (.01202) (.01586) 
7.0 30873 -.18768 -.26080 -.04054 335 
COIS) (.01709) (.01481) (.03720) 
10.0 .52842 90623 -.30765 -.47076 184 
(.02135) (02773) (.01864) (.08035) 
[25 OO ele es oul -.8509 1 160 
(.03012) (.04520) (.02174) (.06960) 
15.0 1.01691 2.06055 -.34404 -.92400 158 
(.04147) (.07748) (.02499) (.05621) 


To study the effects of censoring for this model, a simulation experiment is done 
in which the exponential censoring times has a mean of 1/c=1000. The other 
parameters are p)=1, p,=1, p,;=1, and 8@=0.5 as before. The number of individuals 1s 


N= 50. The results are shown in Table XV. The average relative bias of the parametric 
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estimate is higher for all times than Table VII using the correct model. The average 
relative bias of the nonparametric estimates are about the same as Table VII. Again, 
even with limited censoring, the average relative bias of the renewal estimate has 


computational problems for large ¢. 


TABLE AV 
AVERAGE RELATIVE BIAS 


Hypoexponential Model N=50 (Limited Censoring, c=0.001) (S00 Reps) 

















’ _A a a” A 
Time P(t) P (t) P (t) P(t) #KM 
5 -.07105 -.00384 -.03344 00197 500 
(:00060) (:0008 1) (:00157) (:00136) 
1.0 -.08365 -.03442 -.04223 00167 500 
(00114) (:00156) (00170) (00226) 
2.0 -.06224 -.19090 -.02480 00162 500 
(:00224) (:00279) (:00254) (:00391) , 
5.0 -.00469 -.70151 -.00948 00063 500 
(00561) (:00294) (:00610) (:00789) 
7.0 02631 -.80405 -.01056 -.00350 . 500 
(00798) (00221) (:00856) (01105) 
10.0 07973 -.75567 -.00598 00391 500 
(:01182) (:00230) (01231) (01615) 
12.5 13071 -.59906 00328 -.01794 500 
(:01535) (.00345) (:01553) (:02164) 
15.0 18831 mvle 01761 -.02465 500 
(:01925) (:00598) (:01889) (02797) 


Two additional simulations are done using different hypoexponential distributions 
for the sojourn time in state 1. For these simulations the number of individuals is set at 
N=50 for comparative purposes. The first simulation uses a hypoexponential 
distribution of Uo et p,=0.1, 8=0.5, and c=0.5. Table XVI shows the 


average relative bias and standard error for this model. 
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TABLE xXVvij 
AVERAGE RELATIVE BIAS 


Hypoexponential Model N=50 (p,;=0.1) (500 Reps) 




















A >» A “A 

Time P(t) P (t) P(t) P(t) #KM 

S -.01070 00048 -. 4146] -.Q00004 S500 
(:00036) (:00033) (01923) (:00049) 

1.0 -.01378 -. 00012 -.49394 -. 00087 500 
(:00071) (:00079) (101835) (:00108) 

2.0 -.00851] -.00275 -. 59801 -.00309 500 
(00141) (:00172) (:01730) (:00208) 

5.0 03290 -.00058 -./4519 -.00397 492 
(:00359) (:00541) (:01537) (:00716) 

7.0 .06534 .01743 -./9285 01353 372 
(:003512) (:00817) (:01446) (:01382) 

10.0 11795 08352 -.83604 -.08234 165 
(:00755) (101168) (:01349) (04149) 

[25 .16554 179472 -. 85869 -.39425 85 
(:00971) (:01464) (:01295) (07374) 

15.0 .21680 28205 -.87462 -.69274 60 
(:01203) (01784) (:01258) (:08075) 


Surprisingly, in Table XVI, the average relative bias of the. parametric estimate is 
significantly better than Table XIII. The survivor function of the first passage time to 
state 0 for the semi-Markov model having the sum of two exponentials with mean |] 
and 10 for the sojourn time in state 1 was computed. It was compared to the 
corresponding survivor function of the Markov model of Chapter II having 
exponential sojourn time in state 1 with mean I}. The parameters p,=1, and 0=0.5 
are as before for both models. For large ¢, the two survivor functions are 
approximately P{D>t}=exp[-0.045t] for the semi-Markov model and 
P{D> t} = exp[-0.043t] for the Markov model. Thus it appears that the small average 
relative bias for the parametric estimate in Table XVI is due to the closeness of the 
survivor functions for the two models. The average relative bias of the renewal 
estimate is significantly better than Tables XIII and IV. The average relative bias of 
the asymptotic estimate is significantly worse than Tables XIII and IV. A possible 
explanation for this is that with a mean sojourn time in state | of approximately 11 


and a mean censoring time of 2, the process either jumps to state 0 at first transition or 
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becomes censored due to the expected long sojourn time in state 1. Therefore, the 
Kaplan-Meier estimates for F, will probably contain most of the probability mass at 
small times and relatively little mass for large times due to censoring, causing the 
Kaplan-Meier estimates of F, to be unreliable for large times. The renewal and 
Kaplan-Meier estimates have approximately the same relative bias for tS7. The 
average relative bias of the Kaplan-Meier estimate of the first passage time to state 0 is 
about the same as Tables XIII and IV. For t2/2.5, the Kaplan-Meier estimate is 
significantly better than Tables XIII and IV, however, of the 500 replications only 85 
are still defined by t= /2.5 in Table XVI. 

The next simulation experiment uses a hypoexponential distribution for the 
sojourn time in state 1 that very closely resembles the exponential distribution used in 
Table IV. For this simulation, the parameter are p,=1, p,=1, p,;= 100, §=0.5, and 
c=0.5. Again N=50 for comparison purposes. Table XVII show the average relative 
bias and standard error of the estimates at selected times. The average relative bias of 


all of the estimates are about the same as in Table IV as expected. 


TABLE XVII 
AVERAGE RELATIVE BIAS 


Hypoexponential Model N=50 (p,= 100) (500 Reps) 




















. A A “A A 
Time P Ct) P(t) P(t) P,{t) #KM 
S -.00825 -.02784 28h] -.00435 500 
(.00213) (.00266) (.00574) (.00334) 
1.0 -.01009 -.11168 -.13732 -.00655 500 
(.00387) (.00420) (.00743) (.00473) 
20 -.01009 -.30050 -. 18659 -.00933 500 
(.00684) (.00717) (.01074) (.00880) 
WU) .02098 -.05434 -.31762 02698 333 
(.01543) (.01872) (.01803) (.03958) 
7.0 06513 .52736 -.36331]1 -.18712 213 
(.02182) (.03064) (.02223) (.08288) 
10.0 .16666 1.90635 -.38933 -.61354 166 
(.03337) (.07764) (.02940) (.11659) 
12.5 .28695 4.30938 -.38193 -./9971 161 
(.04577) (.16887) (.03724) (.14178) 
15.0 .44578 699199 -.35141 - 1.00000 160 
(.06199) (.36461) (.04795) (.00000) 
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l) 


2) 


3) 


4) 


5) 


6) 


IV. CONCLUSIONS 


From the results of Chapter III, it can be concluded: 


The maximum likelihood estimate uses the most assumptions about the model. 
It understandably does well when the model used is correct. It is the most 


sensitive to incorrect model assumptions. 


The renewal estimate and asymptotic estimate are biased by censoring of the 
last sojourn time in a state which makes the Kaplan-Meier estimate undefined. 
Further analysis could be done to investigate reasonable methods to make the 


Kaplan-Meier estimate honest. 


The asymptotic estimate has the smallest average relative bias for large times, 
t= 15. However, the bias is always negative. Further analysis could be done to 


find a bias corréction forat 


The Kaplan-Meier estimate of the first passage time to state 0 uses the least 
knowledge. It does well for small times and moderate to large numbers of 
individuals. The Kaplan-Meier estimate and the renewal estimate appear to do 


about as well for small ¢. 


The larger the number of individuals the smaller the average relative bias 1s for 


all the estimates. 


The renewal estimate requires a great deal of computation. In view of the 
simulation results, one recommendation is to use the Kaplan-Meier estimate as 
long as not too many observations are censored, and then use the asymptotic 
estimate for larger times. The asymptotic estimate needs to be used with 


caution if the last sojourn times in state | or 2 are censored. 
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